1. A Quick Introduction to the Lorenz Equations

Chaos is aperiodic long-term behavior in a deterministic system that exhibits sensitive dependence on initial conditions. $\newcommand{\pa}[1]{\left(#1\right)}$ $\newcommand{\br}[1]{\left[#1\right]}$

Any study of chaos begins with the Lorenz equations:

\begin{align*} \dot{x} &= \sigma\pa{y-x}, \\ \dot{y} &= x\pa{r-z} - y, \\ \dot{z} &= xy-bz, \end{align*}

for parameters $\sigma, r, b > 0$: $\sigma$ is the Prandtl number, $r$ is the Rayleigh number, and $b$ is nameless. As Lorenz did, we will study the particular case $\sigma = 10$, $b = 8/3$, and $r = 28$. Note: From here on forward I will be referring to the Lorenz equations as system (1).

In this code, we seek to plot trajectories in three dimensions and create the strange attractor, known as the Lorenz attractor or Lorenz butterfly (which has a fractal dimension between 2 and 3).

2. Simple Properties of the Lorenz Equations

2.1 Fixed Points

Simply from glancing at the equations in system (1) we see the following fixed point: $\pa{x^\ast, y^\ast, z^\ast} = \pa{0,0,0}$, which holds $\forall\,\sigma, r, b$. For $r>1$, another fixed point can be quickly discovered:

\begin{gather*} 0 = \sigma\pa{y-x} \Rightarrow y^\ast = x^\ast; \\ 0 = x\pa{r-z} - y = x\pa{r-z-1} \Rightarrow\boxed{z^\ast = r-1}; \\ 0 = xy-bz = x^2-bz \Rightarrow\boxed{x^\ast = y^\ast = \pm\sqrt{b\pa{r-1}}}. \end{gather*}

We will call these fixed points $C^+$ and $C^-$. As $r\to 1^+$, $C^+$ and $C^-$ converge to the origin to form a pitchfork bifurcation.

2.2 Stability of the Origin

2.2.1 Linear

The linearization at the origin is

\begin{align*} \dot{x} &= \sigma\pa{y-x}, \\ \dot{y} &= rx - y, \\ \dot{z} &= -bz, \end{align*}

obtained by omitting nonlinearities in system (1). Since the third equation is now decoupled, we see that the solution is $z(t) = z_0\exp{\pa{-bt}}$, where $z_0$ is the initial condition. The other degrees of freedom are governed by the system

$$\begin{pmatrix} \dot{x} \\ \dot{y} \end{pmatrix} = \begin{pmatrix} -\sigma & \sigma \\ r & -1 \end{pmatrix}\begin{pmatrix} x \\ y \end{pmatrix},$$

with $\tau = -\pa{\sigma + 1} < 0$ and $\Delta = \sigma\pa{1-r}$. If $r>1$, $\Delta < 0$ $\Rightarrow$ the origin is a saddle point. If $r<1$, the origin is a sink. Note that $\tau^2 - 4\Delta = \pa{\sigma +1}^2 - 4\sigma\pa{1-r} = \pa{\sigma - 1}^2 + 4\sigma r > 0$, so the origin is a stable node for $r<1$.

2.2.2 Global

To show that the origin is globally stable, we construct a Lyapunov function. Consider $V: \mathbb{R}^3\to\mathbb{R}$, defined by

$$V\pa{x, y, z} \equiv \frac{1}{\sigma}x^2+y^2+z^2.$$

Then,

\begin{align*} \dot{V} &= 2\pa{\frac{1}{\sigma}x\dot{x} + y\dot{y} + z\dot{z}}, \\ \frac{1}{2}\dot{V} &= x\pa{y-x} + y\br{x\pa{r-z} - y} + z\pa{xy-bz} \\ &= \pa{1+r}xy - x^2 - y^2 - bz^2 \\ &= -\pa{x-\frac{1+r}{2}y}^2 - \br{1-\pa{\frac{1+r}{2}}^2}y^2 - bz^2, \end{align*}

we see that $\dot{V}$ is strictly negative $\forall\, x,y,z$, so long as $r<1$. $\dot{V} = 0$ only when $\pa{x,y,z} = \pa{0,0,0}$, therefore the origin is globally stable for $r<1$.

2.3 Stability of $C^+$ and $C^-$

Suppose $r>1$ so that $C^+$ and $C^-$ exist. We can determine their stability from the eigenvalues of the Jacobian matrix for the system.

$$J = \begin{pmatrix} \partial_x\dot{x} & \partial_y\dot{x} & \partial_z\dot{x} \\ \partial_x\dot{y} & \partial_y\dot{y} & \partial_z\dot{y} \\ \partial_x\dot{z} & \partial_y\dot{z} & \partial_z\dot{z} \end{pmatrix} = \begin{pmatrix} -\sigma & \sigma & 0 \\ r-z & -1 & -x \\ y & x & -b \end{pmatrix} = \begin{pmatrix} -\sigma & \sigma & 0 \\ 1 & -1 & \mp\sqrt{b\pa{r-1}} \\ \pm\sqrt{b\pa{r-1}} & \pm\sqrt{b\pa{r-1}} & -b \end{pmatrix},$$

after replacing $x, y, z$ with $x^\ast, y^\ast, z^\ast$, respectively. Then,

\begin{align*} 0 &= \text{det}\pa{J-\lambda I_3} \\ &= \text{det}\begin{pmatrix} -\sigma-\lambda & \sigma & 0 \\ 1 & -1-\lambda & \mp\sqrt{b\pa{r-1}} \\ \pm\sqrt{b\pa{r-1}} & \pm\sqrt{b\pa{r-1}} & -b-\lambda \end{pmatrix} \\ &= \lambda^3 + \pa{\sigma + 1 + b}\lambda^2 + b\pa{\sigma + r}\lambda + 2b\sigma\pa{r - 1}, \end{align*}

let $\lambda = i\omega$, so that

\begin{equation} \tag{2} 0 = -i\omega^3 - \pa{\sigma + 1 + b}\omega^2 + b\pa{\sigma + r}i\omega + 2b\sigma\pa{r - 1}. \end{equation}

Both the imaginary and real parts must be zero, we see from the imaginary component that

\begin{align*} 0 &= -i\omega^3 + b\pa{\sigma + r}i\omega, \\ \therefore \omega &= 0, ~ \omega^2 = b\pa{\sigma + r}; \end{align*}

but $\omega \neq 0$ because then Equation (2) would not hold, so $\omega^2 = b\pa{\sigma + r}$ is the only solution. We substitute this result into the real component

\begin{align*} 0 &= - \pa{\sigma + 1 + b}\omega^2 + 2b\sigma\pa{r - 1} \\ &= - \pa{\sigma + 1 + b}\br{b\pa{\sigma + r}} + 2b\sigma\pa{r - 1} \\ &= -\pa{b + b^2 - b\sigma}r - \pa{b\sigma^2 + b^2\sigma + 3b\sigma}, \\ \therefore r &= \sigma\pa{\frac{\sigma + b + 3}{\sigma - 1 - b}}; \end{align*}

note that $r>0 ~ \Rightarrow ~ \sigma > 1 + b$. Let $$r_H \equiv \frac{\sigma\pa{\sigma + b + 3}}{\sigma - 1 - b},$$ where we use the subscript $H$ to denote $C^+$ and $C^-$ losing their stability in a Hopf bifurcation at this value.

3. Coding It Up

# For math operations and variable storing
import numpy as np

# To solve the differential equations
from scipy.integrate import odeint

3.1 Problem Setup

We start by defining all of our variables. In particular, we define $\sigma = 10$, $r = 28$, and $b = 8/3$, as mentioned above. We also define the number of steps in our iteration with the steps variable. Our initial conditions are stored in the X0 and X1 variables, and t is our time array.

sigma = 10
r = 28
b = 8/3

steps = 10000

t = np.linspace(0, 100, steps)

X0 = [0.0, 1.0, 1.05]

# Compute the fixed points
x_fixed = [0, np.sqrt(b * (r - 1)), -np.sqrt(b * (r - 1))]
y_fixed = [0, np.sqrt(b * (r - 1)), -np.sqrt(b * (r - 1))]
z_fixed = [0, r - 1, r - 1]

We now define a function lorenz which describes our system of equations in group (1). It takes in the initial conditions and all parameters in the Lorenz equation.

def lorenz(X, t, sigma, r, b):
    x, y, z = X

    # Lorenz equations
    dxdt = sigma * (y - x)
    dydt = x * (r - z) - y
    dzdt = x * y - b * z
    
    return [dxdt, dydt, dzdt]

3.2 Solving the Differential Equations

To solve the differential equations we use odeint. This decision was made in order to improve accuracy and reduce computational complexity. From the figures that I've generated, I prefer the ones produced by odeint over standard numerical techniques like Forward Euler.

solution0 = odeint(lorenz, X0, t, args=(sigma, r, b))

# Extracting components from the solution
x = solution0[:, 0]
y = solution0[:, 1]
z = solution0[:, 2]

3.3 Plotting

In this section we plot the Lorenz attractor for our parameters. Additional plots can be found in my github.

3.3.1 Using plotly

The following cell imports all of the necessary packages to plot using plotly.

# Custom styling for figures
import sys
sys.path.append('../')
import style

# To create the animations
import plotly.graph_objects as go
import plotly.io as pio
pio.templates["custom_template"] = style.template
pio.templates.default = "custom_template"

This cell sets the automatic perspective to face the front of the Lorenz attractor.

camera = dict(
    eye=dict(x=1.5, y=-1.5, z=0.15)
)

The following cell prints the Lorenz attractor for a single data set.

# Define labels for fixed points
labels_fixed = [
    f'Origin: ({x_fixed[0]}, {y_fixed[0]}, {z_fixed[0]})',
    f'C+: ({x_fixed[1]:.2g}, {y_fixed[1]:.2g}, {z_fixed[1]})',
    f'C-: ({x_fixed[2]:.2g}, {y_fixed[2]:.2g}, {z_fixed[2]})']

fig = go.Figure(data=[
    # Plotting the trajectory
    go.Scatter3d(
        x=x, y=y, z=z,
        mode='lines',
        line=dict(width=0.5),
        name='Trajectory'
    ),
    # Plotting the fixed points
    go.Scatter3d(
        x=x_fixed, y=y_fixed, z=z_fixed,
        mode='markers',
        marker=dict(size=5, color="#9FFFCB"),
        text=labels_fixed,
        hoverinfo='text',
        name='Fixed Points'
    )
])

fig.update_layout(
    title=dict(
        text=f"Lorenz Attractor for sigma = {sigma}, \
r = {r}, and b = {b:.2g}",
        x=0.5,
        xanchor='center'
    ),
    scene=dict(
    xaxis_title='x axis',
    yaxis_title='y axis',
    zaxis_title='z axis',
    camera=camera
    ),
    legend=dict(
        orientation='h',
        yanchor='bottom',
        y=-0.3,
        xanchor='center',
        x=0.5
    ),
    height=700
)

fig.show()  # To print the figure in the notebook

References

Strogatz, Steven. "Lorenz Equations." Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering. CRC Press, 2018. pp. 309 - 328.